//--------------------------------------------
//
//  Test data structure for meshing.
//  THIS IS ONLY A TEST!  DO NOT USE THIS FOR PRODUCTION
//  CODE!!!!
//--------------------------------------------

#include <boost/multi_array.hpp>
#include <CGAL/Bbox_3.h>
#include <set>
#include <vector>
#include <string>
#include <fstream>
#ifndef XDM_DATA_H_
#define XDM_DATA_H_


namespace CGAL
{
  namespace XDM_test
  {
    struct TestPoint
    {
      TestPoint( float x_, float y_,
                 float z_, float w_ ):
        x( x_ ), y( y_ ), z( z_ ), w( w_ ){}
                 
      float x;
      float y;
      float z;
      float w;
    };
    
    class FIndicator : public std::unary_function<int, double>
    {
      const int label;
    public:
      FIndicator(int i) : label(i) {};
      double operator()(int x) const { return (x == label) ? 1. : 0.; }
    };

 static const int nVertices = 8;
      static const int nCubeDimension = 3;
      static const int VertexOrderGrayCode[ nVertices ][ nCubeDimension ] =
        {
          { 0, 0, 0 },  // 0
          { 1, 1, 1 },  // 1
          { 0, 0, 1 },  // 2
          { 1, 1, 0 },  // 3
          { 0, 1, 0 },  // 4
          { 1, 0, 1 },  // 5
          { 0, 1, 1 },  // 6
          { 1, 0, 0 }   // 7
        };
    
    class XDM_Data
    {
    public:
      typedef CGAL::Bbox_3 Bbox_3;
      
    public:
      
      Bbox_3 compute_bounding_box() const
      {
        
        float dI = fLengthI / float( 100 );
        float dJ = fLengthJ / float( 100 );
        float dK = fLengthK / float( 100 );
        
        return Bbox_3( -dI, -dJ, -dK,
                       fLengthI * nPixelI + dI,
                       fLengthJ * nPixelJ + dJ,
                       fLengthK * nPixelK + dK );
       
      }
      
      int trilinear_interpolation( float x, float y, float z ) const
      {
        return TrilinearLabelInterpolation( x,  y, z );
      }
    
      int   operator()( float x, float y, float z ) const
      {
        int nI = static_cast<int>( x / fLengthI );
        int nJ = static_cast<int>( y / fLengthJ );
        int nK = static_cast<int>( z / fLengthK );
        if( InBound( nI, nJ, nK ) )
          return oData[ nI ][ nJ ][ nK ];     
        else
          return 0;
      }

      void Set( int nI, int nJ, int nK, int nValue )
      {
        oData[ nI ][ nJ ][ nK ] = nValue;
      }
      void  Resize( int nI, int nJ, int nK )
      {
        nPixelI = nI;
        nPixelJ = nJ;
        nPixelK = nK;
        oData.resize( boost::extents[nI][nJ][nK] );
      }
      void  SetPixelSize( float fI, float fJ, float fK )
      {
        fLengthI = fI;
        fLengthJ = fJ;
        fLengthK = fK;
      }

      //-----------------------------------------
      //  GetBndPoints
      //-----------------------------------------
      std::vector< TestPoint > GetBndPoints() const
      {
        return oBndPoints_;        
      }
      
      void ReadBndPoints( const std::string & oFilename )
      {
        std::ifstream is;
        is.open( oFilename.c_str() );


        int count = 0;
        int tmp = 0;
        is >> std::skipws
             >> count >> std::skipws
             >> tmp >> std::skipws
             >> tmp >> std::skipws
             >> tmp;
          
        for( int i = 0; i < count; i ++ )
        {
          TestPoint t(0, 0, 0, 0);
          is >> std::skipws
             >> t.x >> std::skipws
             >> t.y >> std::skipws
             >> t.z >> std::skipws
             >> t.w;
          oBndPoints_.push_back( t );
          std::cout << t.x << " " << t.y << " " << t.z << " " << t.w << std::endl;
        }
        std::cout << oBndPoints_.size() << std::endl;
        is.close();
      }
      
    private:
      typedef boost::multi_array< int, 3 > GridType;
      GridType oData;

      float fLengthI;
      float fLengthJ;
      float fLengthK;
      int nPixelI;
      int nPixelJ;
      int nPixelK;
      std::vector< TestPoint > oBndPoints_;
     
      
      bool InBound( int nI, int nJ, int nK) const
      {
        return nI >= 0 && nI < nPixelI
          && nJ >= 0 && nJ < nPixelJ
          && nK >= 0 && nK < nPixelK ;
      }

      
      //
      //  linear interpolate f(x) at x between end points x0 and x1
      //  Note that the interval is [x0, x1), x0 <= x < x1,
      //  x0 < x1
      float Interpolate( float x0, float x1,
                         float f0, float f1, float x) const
      {
        float f = ( x - x0 ) * f1 + ( x1 - x ) * f0  ;
        return f;
      }

      //-------------------------------------------------------------------
      //  interpolate given a label, find the indicator function
      //  value at x, y, z
      //
      //-------------------------------------------------------------------
      float InterpolateLabelIndicator( float x, float y, float z, int nLabel ) const
      {
        double pValues[nVertices];
        float cX = x / fLengthI;
        float cY = y / fLengthJ;
        float cZ = z / fLengthK;
        
        int nX = static_cast<int>( cX );
        int nY = static_cast<int>( cY );
        int nZ = static_cast<int>( cZ );

        FIndicator FLabelIndicatorFn( nLabel );
        
        for( int i = 0; i < nVertices; i ++)
        {
          int nNgbX = VertexOrderGrayCode[i][0] + nX;
          int nNgbY = VertexOrderGrayCode[i][1] + nY;
          int nNgbZ = VertexOrderGrayCode[i][2] + nZ;

          if( InBound( nNgbX, nNgbY, nNgbZ ) )
            pValues[i] = FLabelIndicatorFn( oData[ nNgbX ][ nNgbY ][ nNgbZ ] );
     
        }

   
        // numbering comes from gray code
        //  (look at the definition of VertexOrderGrayCode)
        float fx00 = Interpolate( nX, nX + 1, pValues[0], pValues[7], cX );
        float fx01 = Interpolate( nX, nX + 1, pValues[2], pValues[5], cX );
        float fx10 = Interpolate( nX, nX + 1, pValues[4], pValues[3], cX );
        float fx11 = Interpolate( nX, nX + 1, pValues[6], pValues[1], cX );
        float fxy0 = Interpolate( nY, nY + 1, fx00, fx10, cY );
        float fxy1 = Interpolate( nY, nY + 1, fx01, fx11, cY ); 

        return Interpolate( nZ, nZ + 1, fxy0, fxy1, cZ ); 
      }

      //--------------------------------------------
      //  Trilinear interpolation of labels using
      //   an indicator function
      //--------------------------------------------
      int TrilinearLabelInterpolation( float x, float y, float z ) const
      {

        if( x < 0 || y < 0 || z < 0 )
          return 0;
        
        int nX = static_cast<int>( x / fLengthI );
        int nY = static_cast<int>( y / fLengthJ );
        int nZ = static_cast<int>( z / fLengthK );
        std::set<int> oLabels;
        for( int i = 0; i < nVertices ; i ++ )
        {
          int nNewX = VertexOrderGrayCode[i][0] + nX;
          int nNewY = VertexOrderGrayCode[i][1] + nY;
          int nNewZ = VertexOrderGrayCode[i][2] + nZ;

          if( InBound( nNewX, nNewY, nNewZ ) )
            oLabels.insert( oData[ nNewX ][ nNewY ][ nNewZ ] );
          else
            return 0;
        }
        if( oLabels.size() == 1)
          return *oLabels.begin();
        

        float fBestValue = 0;
        int nBestLabel   = 0;
        for( std::set<int>::const_iterator pIter = oLabels.begin();
             pIter != oLabels.end(); ++ pIter )
        {
          float fCurVal = InterpolateLabelIndicator( x, y, z, *pIter );
          if( fCurVal > fBestValue )
          {
            nBestLabel = *pIter;
            fBestValue = fCurVal;
          }
        }

        return nBestLabel;
      }
      
    };
    
  }
  
}







#endif
